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Image processing mefliod and apparatus 

The present invention relates ta a method and a]^>aratus for processing 3D datasets, 
particularly but exclusively 3D smfacerepresentations. 

A major field of application of fhe invention is the processing of 3D surface 
representations of the himian body, particularly faces. It is often desirable to 
re^ster 3D surfaces having at least a partial overlap eg to combine adjacent 
oveilapping surface regions (eg as obtained by an optical scanner) into one overall 
surface or to compare pre- and post operative images of patients undergoing surgery. 

An algorithm which is conventionally used for registration of 3D surfaces is the 
It^tive Qosest Point (ICP) algorithm, originated by Besl and McKay (P. J. Besl, N. 
D. McKay. A Method for Registration of 3-D Shapes. IEEE Trans Pattern Analysis 
and Machine Intelligence pp.239-256, VoLM» no.2, February 1992 which is 
incorporated herein by reference). Before describing thi^ algorithm the general 
problem of registering 3D surface will be outlined wi/ch refoience to Hgures 1 and 



Figure 1 shows a "^registaf surface SreG which is shown in die form of 
atriangulated mesh of points ql, qj2»...ql3....qm (in general, qi) which is to be rotated 
in three dimendcHis by a rotaj^on matrix R and translated in three dimrasicms by a 
translation vector T so as to merg^ optimally with a '^reference'' surface SreF 

similarly comprising a sinq>lex mesh d points pl» p2» p3 pl3 (in general, pi) and 

having a re^on which ov^aps widi and corresponds to a region of the regist^ 



The proUrai is essentially to find R and T which mimmise the aggr^te distances 
between a set of points on SreG corresponding closeest points on SreF* 

Accordingly the expresdon: 



is minimised irr the ICP algorithm (N being the number of points in the region of 
overiap and qi* being the closest point on SrEG to qi ). Thus the sum of the squares 



surface. 




1 = 1 



(1) 
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of the distances is minimised to find R and T. This is illustrated in Hgore 2 which 

shows die operation of matrix R and vector T on surface SreG» A« resulting 
aggregate of distances di i, d22. ^33 — dii being minimised. In practice the sum of 
the squares of the distances in the region of overlap will never be zero because of 
measurement inaccuracies. 



The efficacy of a standard ICP approach is limited by the following factors: 
IQ a) erroneous data can detrimentally influence the registration; 

b) the algorithm is susceptible to local minima - hence it requires initial alignment to 
be very close and it is hard to consistendy reproduce similar results, and 

c) although die speed and performance of the ICP can be optimised by sanq>Iing a 
15 specific subset of points from the ^register' surface, the optimal selection of points is 

extremely difficult to determine and to lepeaL 

More sophisticated Versions of die ICP algorithm have been devised, eg by Mauier 
et al (Calvin R Nfaurer et al, IEEE Tram Medical Imaging Vol 15 No 6, December 
1996 pp 836 to 848, incorporated herein by reference) who apply a weighting factor 
approximately proportional to the size of the triangles (of the triangulated reference 
surface). Also the registration is based on matching of surfaces as well as points. 
Mauier et al examine the closest vertex on the reference surface and examine eachof 
its adjacent triangles until they determine the triangle to which the vertex projects. 
They claim that this approach produces the actual closest point on the surface 
^ **approximately 993% of the time'* and that "the occasional errors are snail and do 
not affect the accuracy of the registration'*. The nature or cause of the errors is not 
described. 
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Z. Thvaig in "Itmtive Point Matching for Registration of Free-Form Curves and 
30 Surfaces'* International Journal of Computer Vision VoLll, No^i, pp.119^152 
(which is also incorporated herein by refermce) discloses an algorithm in which die 
expression (1) above is modified by introducing a weighting factor which takes into 
account uncertainty in the positions of the points. 

35 It is known that point to triangle matohing is important when pomt to point matching 
is inadequate due to a low sampling density of the rrference surface. We have found 
that an implementation of the ICP algcmthm, using point to trian^e matohing only, 
which neglects the exceptional cases of ridges and peaks leads to unstable behaviour, 
inaccuracy, or lack of convergence. 
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An object of the present invention is to overcome or alleviate the above 
disadvantages of the prior art 

In one aspect the invention provides a method of registering two surfaces with 
corresponding surface regions by associating corresponding points of the surfaces 
wherein a first point on one surface facing a region of the other surface is associated 
with a corresponding point in that region by projecting the first point to a slope 
discontinuity (eg a ridge or peak) of that region. 

In a preferred embodiment, if the first point lies over a ridge» an association between 
Ae first point and a perpendicular extending from the ridge to the first point is 
established. 

In a preferred embodiment, if the first point lies over a peak, an association between 
the first point and the peak is estaHished. 

In another aspect the invention provides a method of iteratively registering two 
spatial regions with corresponding parts by repeatedly: 

i) finding corresponding sets of points in the respective spatial regions, 

ii) rdadvely moving the q>atia] r^ons to minimise an aggregate of the distances 
between the coirespcmding pairs of points, and 

iii) finding new sets of corre^[K>nding points, 

wherein a weighting factor which is depend»t on the relationship between local 
surface nrighbouihoods of candidate pairs of corresponding points is applied during 
the iteration to accd^rate the convngrace or enhance the stability of the iteration. 

In one embodiment die spatial regions are surfaces and the local spatial 
neighbourhoods are smface prop^es. For example the weighting factor can be 
inversely related to the angle between the surface normals of the can<fidate 
correspCTiding points, or could be invCTsely related to the distance between the pairs 
of candidate points relative to die average distaiK:e between pairs of candidate 
corresponding pcnnts. Also the wdghting factor could be dependent upon die 
relationship between respective surface curvatures or textures for example* 

In one embodiment die following function is minimised: 
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N 

2 wi||qi'-(Rqi + T)|| 2 (2) 
i=1 

3 where wi is a weighting factor which is dependent on the relationship between local 
surface neighbouihoods of candidate pairs of corresponding points qi* and qi. 

It wiU be noted that the weighting factor utiUsed by Maurer etalwasH non-relative 
weighting factor, ie dq>end©Qt only on local surface neighbouihood of a point, not 
10 on the relationship between local ndghbomfaoods of points on both surfaces. 

In a related aspect the invention relates to a method and apparatus for registering 
semi-rigid or partially rigid surfaces, particularly but not exclusively the human face. 

j3 It is often useful to align to surface whidi are not identical throughout, but share a 
common ri^d baseline. For example, surface data of the face captured prior to 
maxillo-facial surgery and data capture post-operatively might share the common 
rigid base-line of the forehead, as opposed to other regions of the face which might 
have changed sigaificandy. 

20 

Accordingly in a fuitho* aspect the invention provides a method of registering a 
partially rigid or semi-rigid first surface representation of a subject with a second 
surface representation of that subject wherein a relatively rigid portion of the first 
surface representation is selected and used as input to a registration algorithm so as 
to register with a corresponding portion of that surface representation. For example 
25 the algorithm can be the modified ICP algorithm described above. 

We have found that the relatively rigid portions tend to be registered preferentially 
by the algorithm, particulariy when the algorithm is enhanced by utilising wdghting , 
factors such as those disclosed in accordance with the previously mentioned aspect 
30 of the invention. Hence differeDces due eg to suigeiy can be revealed. 

In one onbodiment the me&od for semi-rigid registration ccnmpiises: 

i) selecting (either manually or automatically) a region of the 'register' surface which 
^3 ^>pear5 to be rigid, 

ii) sampling vertices (possibly randomly) from tiie selected region, and 

iu) feeding^tiie sampled 'register' vertices together with the entire reference surface 
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to an ICP based registration aigorittim (pref^bly an algorithm in accordance with a 
previously-mentioned aspect of the invention). 

The corresponding rigid points of the 'reference' surface will be identified 
automatically as the algorithm converges. 

It is not necessary to select the rigid region(s) manually or indeed on die basis of any 
prior knowledge of the surfaces. Instead the algorithm can be tailored to do this by 
the. use of certain w^ghting factors. In a preferred embodiment points from the 
register surface are sampled automatically according to a random distribution which 
is weighted according to the weights dependent on the similarity of local surface 
neighbourhoods of candidate pairs of corresponding points eg in accordance with the 
second-mentioned aspect of the invention. The quality of selection will improve as 
the algOTitfam converges until in its final stages^ only data from rigid regions of the 
surface are sampled. 

In this and other embodiments in order to ensure that the selected region in the 
'register" surface is indeed rigid, a 3D distance map (between the 2 surfaces) may be 
superimposed over one of the surfaces after the registr^pn process is complete. The 
user may then check whedier the points selected were indeed rigid by their projected 
distance to the oth^ surface. Any errors in selection can then be corrected by 
resampling only points which are particularly close to the reference surface, for 
example. This point resampling can be peformed automatically. This process can 
optionally be performed iteratively to optimise the resulting registration. 



^ Preferred embodiments of the invention are described below by way <^ example 
only with ref^enoe to Hgures 1 to 12 of the accompanying drawings, wherein: 

Hguie 1 is a diagranmiatic representation of the movemrat of a triangulated register 
surface to superimpose a region of overlap on the corresponding region of a 
30 tirangulated rrfermce surface; 

Hgure 2 is a diagrammatic representation of the neariy re^stered surfaces of Figure 
1 showing the neariy minimised distances between pairs of corresponding points; 

35 Hgure 3 is a series of diagrammatic tranverse views of two surfaces as they are 
leistered by an ICP algorithm; 



Figure 4 is a flow diagram of the ICP algoridun; 
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Figure 5 is a sketch perspective view of two surfaces being registered by a method in 
accordance with first aspect of the present invention; 



Hgure 6 shows the search space (in which points projecting to a peak or a ridge may 
^ lie) over a convex pyramidal region of the low» surface of Figure 5; 

Figure 7 shows a modification in accordance with the first aspect of the invention of 
the flow diagram of Figure 4; 

10 Figure 8 is a diagiammatic transverse view of two surfaces being registered on the 
basis of a weighting factor based on the angle between normals* in accordance with 
the second aspect of the invention; 

Figure 9 is a diagrammatic transverse view of two surfaces being registered on the 
15 basis of a weighting factor based on the distance between a pair of candidate 
corresponding points , in accordance with the second aspect of the invention; 

Hgure 10 is a plot of die variation of the weighting factor of Hguie 8 with 0, the 
angle between the nonnals from a pair of candidate corresponding points; 
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30 



35 



Figure 11 is a plot of the variaticm of the weighting factor of Figure 9 with the 
distance d between a pair of candidate corresponding points, and 

Figure 12 is a diagrammatic representaticm of a computer-generated display of 
registered semi-rigid facial representations of a human face in accordance with the 
third aspect of the invration. 

Figures 1 and 2 have already been referred to. 

Referring to Hgure 3aX points pa to pf are shown on a first surface and points qa to 
qe are shown on a second surface. The first surface can be regarded as the reference 
surface and the second can be regarded as the register surface, ie tiie smface which is 
moved to regjister with the reference surface, or vice versa* Each s^ of points is 
shown connected (by hairlines) to the closest points of the other set 

These connections can be found by known methods such as those employing a k-d 
binary tree for example. The k-d binary tree (Henderson 1983, Friedman e( a/ 1977) 
positions the data on a binaiy tree in which each point on the tree is connected to its . 
nearest neighbours. 
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The aggregate of the squares of die distances along the above hairlines is minimised 
in accordance with expression (1) above. 

At least one surface (in this case the register surface) is triangulated and, as shown in 
^ Figure 3b) after £q)plying the rotation and translation determined from exjmssion 1 
to the set of points qa etc, new points qa*, qb* qc' etc are generated by projecting 
normals (shown arrow-headed) to die facets. The distances along these normals are 
the distances from points pa to pf to the other surface and new points qa% qb* qc' 
are the closest points to pa to pf* 

10 

The rotation and translation required to minimise the aggregate of the squares of the 
distances p - q' are found from expression (1) and applied to move the register 
surface containing points qa** qb' ^ towards the reference suiface. 

15 

The process is repealed iterativdy until the aggregate distance (or its mean) is below 
a predetermined threshold. 

This process is summarised in Figure 4, which shows a first st^e 100 involving 
finding the closest points q (on the register surface) to the points p (on the reference 

^ surface), a second stage 200 involving minimising the aggregate of the distance pq 
between the conesponding points tofind R and T of expression (1) and a third stage 
400 involving generating new candidate corresponding points; q* (Figure 3b)) 
following rotation and translation of the register surface under R and T. When the 
mean of squared distances is below a predetermined threshold (intermediate step 

^ 300) the iteration tmninates. 

As noted above, we have found that the above iteration often does not converge 
satisfactorily or even becomes unstable and the reason for tbis appears to lie in the 
conventional methods ci generating new candidate points by projecting normals to a 
30 triangulated surface. 

This probl^ is illustrated in Figure 5. A triangulated regist^ surface is composed 
of points ql to q6 (in general, qO and comprises a convex pyramid (ie pcnnting 
towards the other surface) defined by points q4, qS» q6, q7 and (|S. A refmnce 
35 surface (not shown triangulated) comprises points pi to p6. 

It is assumed that the ICP algorithm has minimised the aggregate of distances plql, 
p2(^,.».. piqi to find the rotation matrix R and translation vector T of expres^on 

(1). 
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In Older to find an improved set of candidate closest points ql\ q2% q3% q4\ q5' q6' 
ete on the register surface the conventional procedure is to project normals (shown 
as arrow-headed lines) to a triangular facet surrounding the preceding dosest point 
^ ql, q2 etc. This works satisfactorily for eg q3* which is generated by flie in^-secticm 
of a normal from p3 to facet qlqSqS. ql% q2' and q4' are also generated correctly by 
this procedure. 

However p6 which lies over a ridge q5q6 cannot be projected normally onto a facet 
10 of the register surface and nor can p5 (which lies over the apex q5 of the pyramid 
referred to above). Accordingly the conventional ICP algorithm utilising point to 
surface matching would fsdl to generate q5* and q6* and this is believed to be the 
reason for the instability of the algorithm in many cases. 

15 In accordance with a feature of the invention, q6' is generated coirecfly by 
projecting a perpendicular to the ridge q5q6 and q5* is generated at the same spot as 
q5 merely by asssodating p5 to the apex q5. 

The geometry of the pyramidal portion bounded by q4, q6, q7 and q8 is illustrated in 
^ Hgure 6. An irregular cross-shaped space defmed by the sets of normals from the 
triangular faces of die pyramid diverges from ridges q4q5, q5q8, q5q7 and q5q6 and 
extends upwardly towards the rrference surface. Any point on the lefiracence surface 
lying within die central region Rl of this space has q5 (the ^x of the pyiamid) as 
its closest point. Any point in one of the four regions R2 of this space has a closest 
^ point defined by the perpendicular from the underiying ridge to that point 

Accordingly the modified algorithm in accordance with the invention searches in 
spaces of the type shown in Figure 6 as shown in Rgure 7. It should be noted that 
that pointe qi surrounded by convex regions made up of five or more triangles will 
have five or more branches in the regions of space corresponding to R2 in Figure 6 

^ and that in general, similar reasoning wai apply to regions of convexity in surfaces 
having other types of simplex mesh» defining their facets - eg hexagonal or other 
regular or irregular polygonal facets. What generates the space of Hgure 6 is the 
convexity of the und^Iying surface - if there is convexity in two orthogonal 
directions then a region sinular to Rl is generated and if there is convexity in one 

35 direction only then a region dnular to R2 is generated. 

This aspect of the invention is not limited to ^irfaces. defined by meshes, but is also 
applicable to surfaces defined by continuous curves or genial surfaces containing 
discontinuities. 



9 



Returning to Figure 7» it will be seen that die modified algorithm searches for qi* in a 
space of type Rl ie in a space defined by biconvexity (step 210) and if no such point 
is found searches in a space of type R2 ie in a space defined by convexity in one 
direction (step 220). If no such point is found then the algorithm proceeds to step 
TXXy in which the aggregate of the distances pq is minimised as in the conventional 
algorithm shown in Hgure 4. If step 210 results in a found point (which will be an 
apex) then the existing point q is identified as an apex and is not bulled in the 
subsequent iteration. If step 220 results in a found point (which will lie over a ridge) 
then the perpendicular from the ridge to the point p is constructed to locate the new 
point on the ridge* 

A further modification to the ICP algorithm will now be desoibed with referem^e to 
Figures 8 to 1 1 in conjunction with expression (2X A wdghting factor wi fw use in 
the expres^on (2) is generated as the following function of the angle 6 (in radians) 
betweCT the normals Npi and Nqi corresponding points pi and qi on die respective 
surfaces: 



20 



25 



30 



35 



wi = exp(«02/a2;52) . exp(-l/4a2) (3) 

and this function is shown in Hgure 10. As a result, during the minimisation of the 
aggregate of the distances b^ween die sets of points pi» qi, those pairs of pcnnts 

which have a similariy ori^ted surrounding region (and are correspondingly likely 
to be true matching points) are driven togeth» preferentially, resulting in a fast^ 
convergence of die algorithm. Preferably 0.1 <; 03. 

Another weig}iting factor wi* can be multiplied or odierwise combined widi with wi 
and is defined as follows with reference to Hgure 9: 

wr = l/(P + d) (4) 
where d is distance and P > 0 . 

This weighting factw enhances the effect of the points on die surfaces which are 
already close togetfaa in driving the algoritfamL towards conv^ence^ 

Tyjpical applications of die above aspects of the invend<»i are: 

a. The afigment of pardally ov^apfnng surfaces of the same object, captured from a 
3D scanning device. 
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b. Fusion of multi-modal 3D surface data, e*g. MRI» CT, Surface scan data» etc. 

Rgure 12 shows a personal computer PC programmed to display a representation of 
a patient's face on a display D, both before and after facial surgery. In order to show 
the change in the patients facial appearance the computer is arranged to run a 
re^stration algorithm which is optimised to register the rigid (eg forehead regions) 
of pre- and post-operative 3D optical scans of the patient's face. 

The apparatus is operated as follows: 

i) a forehead region FH of eg the post-operative surface scan which appears to be 
rigid is selected 

ii) vertices from the selected region FH are sampled, and 

iii) the sampled vertices are fed together with the entire pre-operative surface scan 
(the latter acting as the reference surface) to an ICP based registration algorithm 
(preferably an algorithm in accordance with a previously-mentioned aspect of the 
invention). 

The corresponding rigid points of the 'reference* surface are identified automatically 
as the algorithm converges and the altered areas of the patient's face are displayed 
selectively eg under the control of a mouse M or otiher pointing device. 

In order to ensure that the selected region in the ^register* surface is indeed rigid, a 
3D distance map (betwem the 2 surfaces) may be superimposed over one of the 
surfaces after the registration process is complete. The user may then check whether 
the points selected were indeed rigid by th^r projected distance to the other surface. 
Any errors in selection can then be corrected by resampling only points which are 
particularly close to the reference surface, for example. This point resampling can be 
peformed automatically. This process can optionally be performed iteratively to 
optimise the resulting registration. 

Potential sq[)plications of this aspect of the invration include: 

1. Aligning pre and post op data according to fixed baseline (as described above); 

2. Re^stering surfaces with outliers (in which the oullim are discarded because 
they are not related by the common rigid transformation relating themajority of die 
surface points); 
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3. detracting best possible plane of symmetry from a smface. 
This would involve peifoiming the following steps: 

(a) fit a plane through the axis of approximate symmetry of the surface; 

(b) generate a new surface by reflecting all the vertices of the surface 
across the plane; 

(c) select r^ons on either surface which are reasonaUy symmetrical; 

(d) apply semi-rigid registration technique as described above to align the surface 
(from which symmetrical regions have been identified) with the oth& surface; 

(e) find the mid-point b^ween the s&ts of corresponding vertices (i.e. ni^nal and 
reflected) in the chosen symmetrical re^ons» and 

(f) perform a least-squares plane fit through the set of mid points to produce the final 
plane of symmetry. 

The plane of symm^ry cam be used as a basis for measuring asyirunetry in the 
surface, particulailyin the case of the surface of a human face. One such application 
involves measiuing the impaction of an orbit in a pati^t suffering from 
enophthalmus. The plane of symmietry is used to derive the saggital plane and hence 
the coronal plane (which is almost parallel with the front of the face). The centre of 
each orbit is projected onto the coronal plane to find the necessary <fistances from 
the centre of each orbit 
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Claims 

1. A method of registering two surfaces with corresponding surface regions by 
associating corresponding points of the surfaces wherein a first point on one surface 
f adng a region of the other surface is associated with a corresponding point in that 
region by projecting the first point to a slope discontinuity (eg a ridge or peak) of 
that region. 

2. A method according to claim 1 wh^in the slope discontinuity is a ridge. 

3. A method according to claim 2 wherein an association between the first point and 
a perpendicular extending from the ridge to the first point is established. 

4. A method according to claim 1 wherein the slope disconttnmty is a peak. 

5. A method according to claim 4 wherein an association between the first point and 
the peak is established. 

6. A method of iteratively registering two spatial regions with rorresponding parts 
by repeatedly: 

i) finding corresponding sets of points in the respective spatial regions, 

ii) relatively moving the spatial regions to minimise an aggregate of the distances 
between the corresponding pairs of points, and 

iii) finding new sets of corresponding points, 

wherein a weighting factor which is dependent on the relationship between local 
surface neighbourhoods of candidate pairs of corresponding points is applied during 
the iteration to iaccelerate die conv^gence or enhance the stability of the iteration. 

7. A method according to claim 6 whermn the spatial regions are surfaces and the 
local spatial neighbourhoods are surface properties. 

8. A method according to claim 7 wherein the weighting factor is inversely related to 
the angle between the surface normals of the candidate corresponding points: 

9. A method accor^g to claim 7 wherdn the weighting factor is inversely related 
to the distance between the pairs of candidate corresponding points. 
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10. A method according to any of claims 7 to 9 wherein the weighting factor is 
dependent upon the relationship between respective surface curvatures or textures 
for exanq^le. 

11. A method according to claim 7 wherein the foUowing function is minimised: 
2 wi|qi'-(Rqi + T)[ 2 (2) 

where wi is a weighting factor which is dependent on the relationship between local 
surface ndghbouAoods of candidate pairs of corresponding points qi» and qj. 

^ ""^"^^ ^ preceding claim which comprises ninning a 
modified ICP algorithm. - 

13. A method of registering a partiaUy dgid or semi-rigid first surface representation 
of a subject with a second surface repres«»tation of that subject wher«n a relatively 
ngid portion of the fim surface representation is sdected and used as input to a 
registration algorithm so as to regist.^ with a conesponding portion of that surface 
nqvesentaticm^ 

14. A method as claimed in claim 13 as dependent upon any of claims 1 to 12. 

15. A method as daimed in claim 13 or claim 14 whidi comprises: 

i) selecting a region of a register surface whidi appears to be rigid, 

ii) sampling votices Sxom the selected region, and 

iii) f^g the sampled 'registei* vertices together with the entire reference surface 
to an ICP based registiatioa algorithm. 

16 A method as claimed in claim 13 or claim 14 points from a register surface are 
sampled automaticaUy according to a random distribution which is weighted 
accordmg to the weights dependent on the similarity of local surface ndghbouriioods 
of candidate p^ of conesponding points; 

17. A method of surface registration substantially as described hereinabove with 
reference to any of Hgures 5 to 1 1 or Rgure 12 of the accompanying drawings. 
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